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OTS PRICE 


ABSTRACT 


The thermal convection equations for a thin layer of 
fluid are solved numerically as an initial value probLem. 

The calculations include only those nonlinear terms which 

have the form of an interaction of a fluctuation in the velocity 

and temperature with the mean temperature field. In the present 

t 

calculations, the velocity and temperature fluctuations have one 
horizontal wave number, and satisfy free boundary conditions 
on two conducting horizontal surfaces. 

The computed steady state velocity and temperature 
amplitudes show many of the observed qualitative featu: es. 

In particular, the experimentally observed oundary layering 
of the mean temperature field is correctly * procuced. and, at 
large Rayleigh number, the total heat transport rouna to oe 
proportional to the cube root of the Rayleigh numoer .xovided 
the fluctuating temperature and velocity amplitudes have that 
horizontal wave number which maximizes the total neat transport 
However, the heat transport found here for .tee boundaries is 
about three times the experimental value for tv d boundaries. 
The mean temperature gradient can become r.v- e near the 

boundaries for large Rayleigh numbers ant *ar . horj^^ntal 


scale motions. 


The linear stability of the system is also investigated, 
and it is concluded that the stable solutions for all Rayleigh 
numbers investigated (R 4. 10 6 ) have horizontal wave numbers 
which very nearly maximize the total heat transport. 

The stability study also indicates regions in which two or 
more horizontal wave numbers are required to support convection 


t 


I . Introduction 

This paper describes the results of a numerical investigation 
of the thermal convection equations for a thin layer of fluid 
confined between two plates on which free boundary conditions 
are employed. Our theoretical procedure is to include only those 
nonlinear terms which describe the interaction of the mean 
temperature with velocity and temperature fluctuations. That 
is, those terms responsible for eddy viscosity and eddy 
conductivity effects on the turbulence itself are omitted. The 
above eddy terms (hereafter referred to as fluctuating self- 
interactions) are discarded in a physically consistent manner, 

so that no unrealistic behavior results from their omission. 

£ 

The motivation for this research is to discover quantitatively 
to what extent the turbulent convection problem can be comprehended 
without the fluctuating self-interactions . The system of equations 
obtained by deleting these terms corresponds to closing the 
hierarchy of moment equations at the first nontrivial level by 
discarding third order cumulants. The resulting system of 
equations is complete and involves no empirical parameters. 
Moreover, the gross energetics of the flow are preserved. 

The method of numerical solution consists in integrating the 
Fourier amplitudes of the velocity and temperature fields forward 
in time until the steady state is achieved. This procedure has 
the advantage of assuring the stability of the final state provided 
a sufficient range of initial data is sampled. The present 
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calculations , carried out on the IBM 7090, contain one horizontal wave 

:i i- } 

number and enough vertical wave numbers to ensure the elimination of 

i tj ■ %-i 

truncation errors. In the numerical analysis, we have set the Prandtl 
number equal to unity. However, inspection of the equations which omit 
the fluctuating self-interactions shows that the heat transport is not 
a function of Prandtl number, if the system is steady.^ - 

The calculated velocity and temperature fields show many of the 
qualitative features of the experimentally determined fields. In par- 
ticular, at large Rayleigh number, R, the total heat transport is found 

1/3 

to be proportional to R , provided the velocity and temperature fluc- 
tuation fields have that horizontal wave number which maximizes the heat 
transport. However, the heat transport found here for free-boundary 
conditions is three times the experimental value for rigid-boundary 
conditions. Preliminary numerical studies of the rigid-boundary problem-, 
indicate that for large Rayleigh numbers (R~ 10°) the heat transport is 
about a factor of 2.3 smaller than that for f ree-boundaries and therefore 
approximately 50 percent higher than the experiment. Thus, it appears 
that the boundary conditions are quite important in producing the 
experimental heat transport. The system has two additional failings: 
it turns out that the fluctuating amplitudes are steady in time and the ' 
horizontal plan form of the motions is indeterminate. 

The mean temperature gradient at low R closely resembles the 
experimental temperature profiles. At large Rayleigh numbers ;:j?0 

(R ~ 10^), the gross features of the temperature profiles are 
correctly predicted by our system. The computed mean temperature 

\ ; . , \ s | "\r;' x t r 

gradients are large in a thin layer adjacent to boundary and are 


quite small in the body of the fluid. The gradients near the 
boundary can become negative for motions of large horizontal 
scale, but remain positive for motions of a sufficiently small 
horizontal scale. 

The stability of the steady state solutions against 
infinitesimal perturbation at other horizontal wave numbers is 
also investigated and the regions of instability are delineated; 
These results closely parallel perturbation results at low 
Rayleigh number and support the idea that the most stable solution 
is near the one transporting the most heat flux (Markus and Veronis 
1958) . A finite amplitude stability study, and the associated 
development of a several-horizontal-wave-number system to steady 
state will be the topic of a future investigation. 

The idea that the important features of the convection 
problem are contained in the system which omits the fluctuating 
self-interactions is implicit in the theory of convection advanced 
by Malkus (1954) . In the original formulation of his theory, 

Malkus sought a maximum for the heat transport subject to the 
constraint that the temperature gradient be positive, and that' 
there be a smallest scale of motion participating in the advective 
heat transport. The smallest scale is supposed to be determined 
by the requirement that it be marginally stable in the presence 
of the mean temperature gradient occuring in the fluid. The 
smallest scale so determined furnished a cut off in the .cosine 
representation of the mean temperature gradient. The assumption 
that the heat transport was maximum under the above constraints 
then led to an explicit form for the temperature gradient. 
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A more recent formulation of the Malkus theory by Spiegel (1962) 
replaces the cosine representation of the temperature gradient 
by an expansion in terms of the set of eigenfunctions, which 
are marginally stable on the mean temperature gradient . This 
version of the Malkus theory is exactly equivalent to the 
system considered here, provided the horizontal scale of the 
motions is such that the mean temperature gradient is everywhere 
positive. In this sense, our numerical results contain, as a " 
special case, the exact solutions to the Malkus theory for one 
horizontal wave number and free boundaries. 

In this connection, a comparison of our computed temperature 
gradients with the predictions of the Malkus theory is relevant. 

In making this comparison, we must keep in mind that the system 
considered here is explicitly confined to only one horizontal 
wave number, whereas Malkus makes ro explicit references to the 

« 

nature of the horizontal-wave-number spectrum. We do not confirm 
_2 

the z law for the gradient outside the boundary layer as 
predicted by Malkus, nor do we find a sharp cut off in the cosine 
spectrum of the temperature gradients. 

At low Rayleigh numbers (R< 2000) our numerical results 
are in agreement with the calculations of Malkus and Veroni:* (1958) 
and Kuo (1961), who have obtained perturbation solutions to the 
connection equations. A procedure similar to ours has been used 
by Salt zman (1961) in studying the complete convection equations 
for R & 6000. Our approach differs from his in that we 
are able to allow very many vertical modes to be excited. 


-5- 


whereas his results are limited to one vertical mode. Our 
results indicate that it is essential to allow many more vertical 
modes than horizontal modes, if large truncation errors are to 
"be avoided. Thus, at R * 4000, 10 vertical modes must be 
included to obtain realistic temperature profiles. 

II . Theory 

a) The Equation of Motion and Boundary Conditions 
We consider a thin layer of fluid confined between two 
infinitely conducting plates located at z = 0 and z = d. The 
lower plate is maintained at zero degrees, and the top plate 
at a temperature -t q# on an arbitrary temperature scale. The 

A 

direction of gravity is specified by the unit vector -k . The 

equations appropriate for our system are the Boussinesq 
2 

approximations to the Navier-Stokes equations. We shall write 
these equations in a form in which the velocity and temperature 
(v and T) as well as the coordinate and time (r and t) are 
nondimens iona 1 . The only parameters entering the equations 

JT 

are then the Rayleigh number, R, and the Prandtl number, o. 

The equations are 

7*v =0 (1) 

'N 

(i ) 7 2 v = — 7x(7x(v-7v) + R7x(7xkT) (2) 

VfT at . y a 


(| - ) T = - 7-(vT) 


(3) 
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Equation (2) is actually the double curl of the momentum equation, 
and lienee the pressure variable is absent. The nondimens ional 
variables are related to the dimensional ones (denoted by primes) 
in the following way: 




T = T'/d 


9 = r'/d 





Here h is the thermometric diffusivity of the fluid. 

The boundary conditions on the velocity field are derived from 
the requirement that the fluid exert no shear on the confining plates 
This, together with the continuity equation, implies that all even 
derivatives of the vertical velocity, w, vanish on the boundary. In 
the nondimens ional notation the boundary conditions are, , 


am 

3Z 


(4) 


- and 

T (0 „ t) = 0; T(l,t) = -l . (5) 

b) Discard of the Fluctuating Self-Interactions 
It is convenient to resolve the temperature field into a 
horizontal mean plus a fluctuating part; 
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T = -z + + ,©(r , t) . (6) 

Here, iji(z,t) is a horizontally averaged distortion or the conduction 
state and 0(r,t) is the fluctuation of the temperature from its 
distorted value. In view of the boundary condition (5), and the 
interpretation of 0 as a fluctuation from the horizontal mean, we. 
may write . 

♦ CO, t) = Ulyt) = 0 (7; 

M(x,y,0,t) = 0(x,y , 1 , t) = 0 = 0 . (8) 

The bar on equation (8) indicates an average over the horizontal. 

We now introduce equation (6) into equations (1), (2) and (3) and 
■ subtract from each of the resulting equations their respective 

I * 

horizontal mean. We find 

(Jj- ^ V 2 w = RV 1 3 @ + i |vxvx(v- Vv)| z 


- v 2 J © = (^1 - ^ ) w - v • (v® - kw©) 

tie -ip;* 


There are two more equations, for the x and y components of the 
velocity field, but these are not necessary for our problem. The 
last terms in the equations above for w and © have the form of a 
deviation .of a bilinear fluctuating quantity from its horizontal mean 
( f luctuating self -interaction) . By discarding these terms we obtain 
the system to be investigated; 
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fi 4 - V3 ) V 3 W = RV 20 (9) 

vcr y 1 

lit - '”0 ® = 3 w (10) 

1 

_ 31 ") ^ = . 3_ , (11) 

^z s y v 

where: s(z) = 1 - ff = - ~ T 

l 

The significance of omitting fluctuating self —interaction can 
be expressed formally by examining the hierarchy of moment equations 
obtained from equations (l)-(3). By multiplying equations (2) and 
(3) by v(t') and T(t') and ensemble -averaging the appropriate sums 
of the resulting equations, we obtain the time evolution equations 
for the correlation coefficients <vnVj'>, and <v j T'>, and <T T'>. 

These equations couple the above second order moments to the transfer 

terms, which are cubic in v and T. 

Since the system contains a non-vanishing first-order moment, 

, the transfer terms contain both correlated third-order moments 

(cumulants) and products of first order moments with second-order 

moments. The discarding of the fluctuating self-interaction then 

corresponds to closing the system of moment equations by discaiding 

3 

the third order cumulants* In the absence of mean fields this 
procedure would be empty. 

We must now verify that our procedure of deleting third-order 
cumulants does not lead to physically unrealistic results. For our 
procedure to be acceptable, the system of equations (9), (10) and (11) 
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must obey the conservation laws associated with the complete set 
of convection equations, and they must be free from unphysical conse- 
quences of . the sort recently discussed by Ogura (1962), for a similar 
problem in isotropic turbulence. Ogura has demonstrated that the 
assumption of zero forth -order cumulants and nonzero third -order 
cumulants is incompatible with a positive energy spectrum for all 
wave numbers. 4 

With regard to the last point, it should be noted that the 
positive definite character of the kinetic energy wave number spectrum 
and the spectrum for the square of the temperature field follows * 
directly from the fact that it is possible to write the equations 
which delete third order cumulants in terms of amplitudes rather than 
moments. We observe that the amplitude equations (9), (10) and (LI) 
all have real coefficients; hence, the square of any amplitude will 
remain positive for all time if the amplitude is initially a real 
number . 

The conservation of entropy and kinetic energy are also 

r 

preserved without the fluctuating self- interactions . By multiplying 
equation (10) by 0, equation (11) by , and adding, we obtain after 
integrating over the entire volume of the fluid , the equation of 
conservation of entropy, 

r f 

\ ft i l 0 l + I* l 3 j + \ l?*| 2 + I v ©| n = {w ©} . 0,2) 

V V 
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Here the v subscript indicates an integration over the entire volume 
of the system. 

We observe that equation (12), with a corresponding one for 
the conservation of the kinetic energy 5 of the flow are exactly 
the same as those with the fluctuating self- interaction included. 
Contributions from the latter may be reduced to surface integrals 
which vanish. 

C) Fourier Decomposition of the Equations 

It is convenient to work with the Fourier components of the 
equations (9) -(11) rather than their space-variable form. The 
free boundary conditions make the sine series appropriate. We 
therefore write : 

w(r,t) = J f a (x,y) w n a sin n 71 z 

n,a 

®(r,t) = £ f a C*>y) 0 n ° sin n 71 z 
n,a 

t(z,t) = sin n n z • 

n 

Here f (x,y) is ah arbitrary set of orthonormal functions generated 
by the operator 7j a , and obeying appropriate periodic boundary 
conditions in the horizontal: 


7 1 3 f a (x,y) = - n 2 a 2 f a (x J y) 
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and 


If 


a f a' ,2 " 6 aa' 


Introducing the above representation into equations (9), (10), 
and (11) gives the following set of equations for the amplitudes 

w n> V “ d V 


(— + n 2 + a 3> ) w a = } ClZ — 3 Ft a 

Va St J n n 2 + a 3 n 


(13) 


(It + + “*) ®n° = “n“ - i I P*p (‘“n+p + " (n -P> “ a |n-p|) (14) 


P = 1 


(It + nS J *n * ¥ I I "pVntp + o( P- n) ^Ip-pl) < l5) 


p=l a 


where 


X = R/tt 4 , 
t = rr 2 t 


and 

ct(x) = 1, x > 0 


Jf 


= 0,, x = 0 
=-l , x < 0 
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Manipulation of the convolution terms in (14) and (15) is aided 
by the following identities: 

l V B n+p + ^"-P^ln-pl 5 = l B p (A |n-p| " A nV ' 

P P 

and, 

l y B m-P + °<P- n)B rp-nl ) “ 1 V Vp + ”<P- n)A lp-nl ) ' • 

p p 


There are two conservation equations derivable from (14) and 
(15). The first is the Fourier representation of equation (12) for 
conservation of entropy. The other is the equation that partitions 
the total heat flux between conduction and convection; and it is 
derived by multiplying equation (15)' by 1/n and summing over n. 


We find 


l (& It * 0 % ■ f I 


(JL) 


a 


n 




a 


n 


( 16 ) 


n 


n,a 


where 


n 


= - TT 


n Y n 


Here the 8 's are the cosine transform of the mean temperature 
gradient. In the statistically steady state, equation (16) is the 
equation for the total heat flux, which is a constant of motion for 
the system. We now introduce a quantity N(t), the total heat flux 


at the lower' boundary: 
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N(t) = 1 + £ B n (t) 
l 


I 

If the mean amplitudes are constant, 


M _ i , n 2 V a n a 
N " 1 + T L ^n ®n 
n,a 


(17) 


V* 

JF 


In our units, the conduction state transports unit heat flux and 
this equation is the nondimensional form for the familiar equation 
for the total heat flux. 

d . Structure of the Equation 

Before proceeding to the numerical results, we give a brief 
resume of the pertinent qualitative features- of the system defined 
by equations (13), (14) and (15). First of all, we note that the 
horizontal wave numbers, a, are coupled only in their effect on the 
mean temperature field \|r. This interaction occurs diagonally in the 
sense that each a "interacts only with itself. As a consequence 
there is a degeneracy in the horizontal plan form of the motion; 
the system is insensitive to the particular cell shape. Moreovei , 
the number of a's is also indeterminate. The simplest situation 
is to have a single a support the motion and we investigate only 
this case here. 

A single a will give nontrivial answers for the amplitudes 
w and 0 only if it lies within a certain range. The range of o 
which will not support convection is obtained by assuming w and © 
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to be small, and demanding that they subsequently decay. If 
w and 0 are small, i|r will be small to second order and our question 
is equivalent to that of marginal stability. The system then will 
not support convection if 


Cl + qg) 3 , R 


a* 


* T7 4 


(18) 


Conversely, we assume that the steady state values of w and 0 will 
be nonzero if a lies in the range complementary to (18). 

The time behavior of the system is complicated by nonlinear 
effects. In the approach to the steady state, our numerical results 
indicate that the system executes overdamped oscillations with an 

i 

ever increasing period of oscillation. This last remark is 

understandable since w and 0 become marginally stable as t - <*>. 

If the mean field, \Jf , is statistically steady as t «>, we 

7 

m £'/ use a theorem of Spiegel to show that w and 0 are independent 
of time. Spiegel has shown that the principle of exchange of stability 
is valid (for free boundaries) in the presence of the mean gradient 
corresponding to the steady state solution to the mean temperature field 
given by (15). This implies that the growth rates for the appropriate 
eigen-function expansion for w and © must all be zero in order for 
there to be a statistically steady state. 


Ill Numerical Procedure 

In performing the numerical integration, we discard from the 
onset those Fourier amplitudes which will be zero in the steady state. 
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We assume that the steady state amplitudes m, 6, and e ’nave even 
parity about the mid point z = 1/2. This means that the even sine 
modes of m and 0, and the odd cosine moaes of e will have zero 
amplitude in the steady state. We therefore put their initial 
values equal to zero. The equations of motion (13), (14) and (15) 
then imply that the odd parity modes will remain zero for all 
subsequent time. Defining "5^ = we ma ^ rewrite (14) and (15) 

in a more convenient form: 

03 

(|_ + n 3 + o 3 j %=%+!£ e p («> n+ 2p + CT( n - 2 p)i»| n . 2 pl> O'") 

P =1 


(§7 + 4n3 ) p n = ’ 2v *™ 2 1 u ’p (0 2n-p + 0 (P- 2n > p '\2n-p\ ) . (15,) 

. p=l 

Here, we have dropped the a superscript since we are interested 
in the system containing only one a. Equation (13) remains 
unchanged and the total heat flux is computed from equation (l7 ) . 

Our procedure for integrating these equations is to assign 
an initial set of amplitudes to uo^, © n> P n , and allow the system 
to evolve to the steady state. In doing so, we must truncate 
the infinite set of equations. Our procedure in this matter is 
to set all amplitudes i» n , © n , 7 of index greater than a certain 
integer, n Q , equal to zero. This method of truncation guarantees 
exact conservation of heat flux and entropy for the abbreviated 
system. Since »\ is generally large and the © n 's decrease rather 
slowly, we see from (15') that amplitudes for P n above n 0 /2 will 
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have significant truncation errors. Truncation errors are assumed 
to he negligible if increasing n 0 does not appreciably alter the 
value of the total heat transport The total ftumber of "g N modes 
iincluded in these calculations ranged from 20 modes at R = 4000 to 
80 modes at R 10 6 . The errors in the total heat transport due to 

4 

the above are estimated to be less than one part in 10 3 

The integration forward in time was continued until constancy 
of heat flux (16) and entropy (12) was achieved to one part in 10 4 . 

The time, in t units, necessary to achieve this ran from ~ 1.4 at 
R = 4xl0 3 to 0.3 at R = 10 s . At high Rayleigh numbers, this 
criterion was not too satisfactory, since constancy of heat flux and 
entropy were achieved long before the amplitudes w and P 1 became 
steady. For these cases, it was necessary to check the time derivatives 
of the slowest evolving amplitudes, w 1 and . The system was 
observed to be steady if the derivative of w 1 was less than 2 percent 

of Wj . 

Examples of the time evolution for the total heat transport 
N(t) are given in Figure 1, for R - 4x1 Of , 10 4 , 10 5 and a - 1-5. 

The system was started in the conduction state at t = 0, with all 
fluctuating amplitudes uj N and 9 N equal to zero, except uj 1 , which 
had an initial value of unity. The convection is seen to develop 
initially by way of large oscillations, and to decay to the steady 
state with overdamped oscillations, whose period becomes increasingly 
larger. The time scale of the initial oscillations in these curves 
is of the order of the growth rate time • in the conduction state. 


-17- 


IV Discussion of Results 

The computed steady state amplitudes are shown in figures 2—13,. 
The normalization for B, w, and ® are given in the captions, while 
T(z) requires no normalization. The graphs of T(z) are in a 
reflected coordinate system to conform with ah accepted procedure. 

The values of a in figures 2-9 were chosen so that the heat 
transport is very near its maximum. We now discuss in some detail 
the physical features of the steady state amplitudes ^ ® r " > end T. 

a ) The Mean Temperature T, and Mean Gradient fl (z ) 

The mean fields, T and 6 in figures 3, 5, 7, 9, 11 and 13 
have an interesting behavior near the boundaries. At low Rayleigh 
number, these fields closely resemble the perturbation results of 
Malkus and Veronis (1958) , but the temperature gradient is slightly 

negative in the central region. As the Rayleigh number is increased, 
the negative temperature region collects closer to the boundary 
while in the central region, the temperature gradient becomes 
extremely small but positive. 

The negative temperature gradient boundary region is apparently 
produced by an overshoot pheonomenon. These occur typically for 
motions of large horizontal scale (small &) and disappear for 
motion of small horizontal scale. (See figures 9, 11 and 13.) If 
the motion has a large horizontal scale, an element of fluid close 
to the lower plate moves in a region of high temperature for a 

When it eventually, turns upward, it moves 


relatively long time. 
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unchecked by eddy processes and penetrates the body of the fluid 
with an excessive heat flux. The mean gradient accommodates this 
motion by turning negative. The negative P region then checks 
the velocity field, so that the advective heat transport decreases 
toward the middle of the fluid. We see here evidence for the 
non-local property of the flow; if the Rayleigh criteria for 
convection were applicable locally, a negative P region would not 
persist in the steady state. 

For motions of small h- rizonta! ._aie (figures 11 and 13 ) 
the situation is somewhat different. In this case, an element 
absorbs little heat from th«. lower boundary region and loses it 
quickly by conduction because it belongs to a vertically elongated 
cell pattern. It also loses momentum by viscous drag, and attains 
its terminal velocity before reaching the central region of the 
fluid (see figures 10 and 12) . To maintain constancy of heat flux 
the central region must conduct rather strongly, so that the 
mean gradient becomes large there. 

b) Velocity and Temperature Fluctuations ^ 

The velocity and temperature fluctuation fields are shown 
in figures 2, 4, 6, 8, 10 and 12. We observe that the velocity 
fields, for all Rayleigh numbers, have an extremely large first 
mode. For example, at R = 4000, (figure 2 ) Wj represents 997 c of 
the total velocity amplitude, while at R = 10 6 (figure 8) w. is 
957o of the total. On the other hand, the © n modes decrease rather 
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slowly as n increases. 

The above behavior of the w n and ® n spectra displays the 

character of the nonlinear coupling in our system. Thus, the Lei.m 

U, e tend to be the dominant contributor to 3 n [see equation (IV)] 

J* It ^ 

for reasonably small n. Conversely,, terms of the form 8 n m-^ and 8 n— 1 ^ 3 . 
tend to be the dominant contributors to © n [equation (14')] v The 
nonlinear coupling scheme in the equations of motion is therefore 
highly nondiagonal, as opposed to the case of isotropic turbulence. 

The strong nondiagonal coupling in the system of Fourier modes 
is a result of the distortion of the mean temperature profile com- 
bined with the pressure and dissipative forces for incompressible 
flow. The above forces are directly responsible for the occurrence 
of sixth-order derivatives in the marginal stability problem, of 
which the steady state amplitudes u) and © are solutions in the 
presence of the mean field 3 - If we solve for the velocit; V odes 
o) n , in the presence of the mean gradient 3, hy using the iteration 
technique of Section V, we see that the higher modes of w n are 
suppressed by a factor ~n -6 . For a reasonable 3/ this factor 
results in the higher uj n modes making only a small contribution to w. 
c ) T emperature Gradient Spectrum 

The cosine spectrum of the mean temperature gradient, 3, - s 
given in figure 14 for R = 10 4 , 10 5 , 10 6 and a = 1.5. We have 
connected the points with a smooth curve for the sake of clarity. 
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We notice a tendency for the lower modes to saturate at R n = 2 , 
which corresponds to the small gradient outside the boundary layer. 

In fact, if = 2 for all n, £(z) is a 6 function, and the gradient 
vanishes everywhere except at the boundary, where it becomes singular 
At large Rayleigh numbers, the p^ spectrum- is nearly Gaussian 
for small n, but decreases more rapidly at large n. 

The tendency for the B n 's (for small n) to approach 2 as an 
'upper bound is closely connected with the fact that the velocity 
field is marginally stable on the mean temperature gradient, p. 

This feature is brought out more clearly by examining the relation 
connecting the mean gradient field, p^ and the Rayleigh number R. 
Using the iteration method of Section V, we find 


a 2 ” R 



_!Z4 X a + a 8 ) 5 , a 

1 - 1 Pj L C(2n+1) 3 + a 2 l 3 ^ e n 

2 1 


®n+l ): ” 


This series for R" ^ converges rather rapidly for all the p 
which have been computedjand the terms explicitly written in 
equation (19) give R to an accuracy of ~ 207, at R = 10 6 . We 
note that for this equation to balance at large R, P a must approach 
2, and the remaining lower modes must decrease rather slowly as 
n increases. 

The computed spectra (figure 14) are qualitatively quite 

8 

different from the one derived by Ma _ku& - '■ His spectrum is given by 


V = 2 ("l - 2H~ri) 


(19) 
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Here 2^+1 is a cut-off in the P n spectrum, an: 5 it is the total 
heat flux in our units. Marginal stability is achieved at a much 
lower Rayleigh number for this spectrum than for the ones computed 
here . 

With regard to the kaikus theory, figures 10 and 11 are relevant. 
For this case (R = 10^, ct. = 6.0) the temperature gradient is „ 
everywhere positive except near the boundaries where it approaches 
zero. The fields in figures 10 and 11 therefore fulfill all the 
requirements of the kaikus theory as formulated by Spiegel (1902) . We 
note for this case that the total heat flux is ~ 22, whereas kaikus 
obtains a heat transport of ~ 11 for free boundaries. In making 
this comparison, one should remember that these computations were 
made for a s.ingle horizontal wave number, whereas the Kaikus 
presumably allowed for a full spectrum of o's. However, if we 
interpret the computed heat transport as an upper bound to the 
heat transport as the kaikus theory prescribes, we conclude that 
for free boundaries the actual upper bound is at least a factor of 
two larger than that obtained by kaikus. 

d ) The Total Heat Transport as a Function of a and R 

The total heat transport, as a function of R and a is given 
ii^ figure 15. The Rayleigh numbers are indicated in the figure. 

These curves closely resemble the perturbation calculations at 
small R, but become increasingly broadened as the Rayleigh number 
increases. For a given R, the heat transport is entirely conductive 


“ fc. 4* 


(N « 1) if a lies outside the bounds prescribed by equation (1?). 

i 

The value of a, which maximizes the heat transport is l/f"2 at the 
critical Rayleigh number (F = 657), and apparently increases 
linearly in R 1/3 , at large R. It is well represented at high 
Rayleigh numbers by the formula: 

1/3 • 

a ^ 0.7 + 0.01 R 
max — 

The data on this point is not entirely conclusive because of the 

large breadth of the curves. It should be pointed out that 

equation (20) cannot be a correct asymptotic formula since a |nax 

is proportional to R 1 ^ 3 and the value of a, beyond which a single, a 

cannot support convection is proportional to R (equation 18). An 

estimate of the Rayleigh number beyond which (20) is incorrect is 

not warranted by the accuracy of the curves, but according to 

1 8 

equations (18) and (20), it is R < 10 

The maximum heat transport as a function of R is given in 

figure 16. For R £ 3000 the data is accurately represented by 

1 /3 

the following R law: 

N = 0.31 R 1/3 (21) 

Experimentally, the Nusselt number N is 0.085 R 1 / 3 , for large 
R, and rigid boundaries (Jakob, 1959). We see no evidence for an 
intermediate R 1 / 4 law, but such a law may only be obtained in the 
rigid boundary problem. Below R -'lO 3 , the data fits smoothly Lo the 
perturbation calculation of Malleus and Veronis (1958) . 
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The discrepancy between (21) and the experiment is parti; 
a result of eddy processes which our procedure oir.its. The use of 
rigid boundary conditions may improve the agreement, but if (21) is 
corrected for boundary effects as done by Malkus (1960) by decreasing 

N by tn, - re remains a discrepancy of a factor of 2. 

If we choose horizontal wave numbers such that the mean gradient 
is everywhere positive (figures 10 and 11) the discrepancy is 
reduced to 1.8. The latter fields, however, have the unattractive 
feature of having a large temperature gradient in the central 
region of the fluid. 

V Linear Stability of t he Fiel Is 

The velocity and temperature fields we have so far discussed 
are stable against the introduction of a disturbance of the i iflL S L 
horizontal wave number for which the fields were computed. Ihis 
stability is inherent in the method of integrating the equations. 

The stability of the steady state amplitudes against disturbances 
at wave numbers a' other than that a which supports the convection 
process has not yet. been assured in our calculations. The question 
of stability of the - solutions against disturbances of finite 
amplitudes leads directly back to the multi-a system of 
equations (13), (14) and (15). We should assume a whole spectrum: 
of q,'s are initially excited, let them evolve to steady state, and 
repeat the calculation for an ensemble of initial conditions. We 
shall be content here with an investigation of the linear stability 
of tlie system. This problem has some intrinsic interest, but our 



main purpose is to lay the framework for an investigation of the. 
rnulci-a system. 

It is convenient to pose the linear stability problem in 
terms of the Fourier amplitudes (equations (13), (14) and (15)). 

We suppose that the system w n a i , > and B n have their steady 

state values, introduce disturbances 5w n a , 6© n and f>8 n and ask 
whether the latter grows or decays. Since and a are not coupled, 
$r must decay initially. The problem then is reduced to determinin 
the growth rates for 6w n a and 6® n a in the presence of the mean 
gradient r^. Since exchange of stabilities has been proved foi 
this system 10 , we know that the system w n a i , © n a i , and e n will be 
stable if the smallest critical Rayleigh number R c for the 
perturbation system, 6 u}a n , 6®°^ i- s larger than the Rayleigh number 
for which w a i and © a i were computed. 

The marginally stable amplitude 6w n , and 6© n satisfy 
equations (13) and (15) at a Rayleigh number R c , with the time 
derivatives put equal to zero. Since the smallest R c takes an 
eigen function even about z = 1/2, we may abbreviate the 
perturbation system by eliminating the even sine modes from the 
velocity and temperature fluctuations. Defining 

.% “ 5w 2n-l 
^n - ®2n 5 

we may eliminate $© n by using the steady state form of equation (13) 
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and write the marginally stability problem in the following matrix form: 


A(P) = u<P 


( 22 ) 


where 


A ( p ) " T ( 2n-l) a +a 3_ i 3 { 6 nm + 2 ^ P, n-m' _ P n+m-l' ) } 


u R » 


In siting the matrix A, we have used the alternative form for the 
convolution term in equation (14). 

The largest eigen value, ji (smallest R c ) , of equation (22) may 
be obtained by the matrix iteration technique (Hildebrand, 1952) . Since 
the first sine mode of the velocity will be largest, we may conveniently 
begin the iteration on a vector containing only this mode. ..Del ating 


ll> = (1, 0, 0, ..., 0 , . . . . ) 


we may write 

S±-. = u = lin (23) 

V“ ) maX n-“ <llA n ‘ 1 |l> 

The convergence of the iteration scheme is quite rapid because of 
the structure of the A matrix. At the highest Rayleigh number 
considered, R = 10 6 , the 11 th iteration gives R c to one part in L0 6 . 
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5 

The calculated R c (rr)'s are shown in figure 17 for R = 5x10 . 

The particular a. which supports the mean field, p, labels the 

3 

various curves. At low Rayleigh numbers, (R < 5.10 ) these curves 
closely resemble those produced by perturbation calculations. 

Above R ~ 10^ they become increasingly distorted; steady state 
amplitudes of small horizontal wave, numbers are enormously unstable 
with respect to an introduction of a disturbance at large a. 

In figures (18) and (19) we give the zones of instability 

4 5 

for the computed amplitudes for R - 10 and R = 5x10 . In these 
graphs, a.} is the wave number that supports the convective process;, 
and o s is the wave number of the perturbation amplitudes. The 
regions of instability are indicated by the shaded areas, whose 
outer boundaries are lines of marginal stability. The line a n = a ? 
is a trivial case of marginal stability. The value of a at which 
the two curves cross represents a solution which is inf initessimally 
stable against all other a's. This value of a begins at 1/V2 
at the critical Rayleigh number and increases slowly with increasing 
R. The rate of increase is seen to be slower than that a which 
maximizes the total heat transport. Referring to figure (15) we 
see that the use of the most stable a instead of will not 

appreciably change the total heat transport. 

The zones which linear stability theory predict must have two 
or more ex's supporting convection are indicated by the cross-hatched 
regions in' figures (18) and (19). These regions are obtained by 
perturbing the a x fields at a s , assuming that the a 2 field 


subsequently dominates the convection, and then demanding that 
the 02 fields be unstable with respect to a perturbation at 
The cross-hatched region is then bounded by the descending mar- 
ginally stable curve and its reflection about the 45° line. At 

3 

small Rayleigh number, R < 4x10 , this area vanishes but it 
gradually increases with Rayleigh number. 

Concluding Remarks 

The temperature and velocity fields computed here with 
the fluctuating self-interactions absent show qualitatively a 
reasonable behavior. The boundary-layering of the temperature 
field, which is found experimentally, is faithfully reproduced by 
the system, and the heat transport has the experimentally determined 
dependence on Rayleigh number. In this respect, our results for the 
velocity and temperature amplitudes, as well as the stability analysis 
of the fields, confirm the original ideas of Malkus . However, our 
result for the heat transport for free boundary conditions does not 
agree quantitatively with Malkus. 

The only disquieting features of the results are the 
negative temperature gradients which can occur near the boundary for 
small cr , and the rather large amount of heat transported by the 
system. • Aside from eddy processes, there are two other modifications 
in tne system which must be explored before its quantitative accuracy 
can oe properly accessed. 

First, the use of the more realistic rigid boundary conditions 
will enable one to examine quantitatively the role of the eo.p processes 
in producing the experimental temperature profile and the total boat 
flux. The presence of shear forces at the boundary will decrease the 
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computec heat flux, and in checking the development of large scale 
horizontal motions there, it will reduce the negative temperature 
gradient. Preliminary indications are that the use of rigid- 
boundary conditions decrease the total heat transport by a 
factor of 2.3. Secondly, the introduction of several horizontal 
wave numbers will make the system more realistic, particularly 
at large Rayleigh numbers. It will also permit a study of finite 
amplitude stability of the system. The above modifications are 
currently under investigation and will be reported on in the near- 


future . 
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FOOTNOT3S 


1. For a discussion of the Prandtl number dependence of the heat 
transport for the complete system see Kraichnan, R. H., 1962: 
Turbulent Thermal Convection at Arbitrary Prandtl Number. 

P hysics of Fluids _5 1374-1389. 

2. See e.g.. Chandrasekhar, S. , 1961: Hydrodynamic and Hydro- 

magnetic Stability. Oxford at the Clarendon Press, p. 16. 

3. Discarding third order cumulants is quite different from dis- 
carding third order moments . The latter procedure has as a 
consequence that no steady state nontrivial amplitudes exist. 

For an investigation of the dynamics of decay for zero third- 
order moments see Deissler, R. G., 1962: Turbulence in the 

. Presence of a Vertical Body Force and Temperature Gradient, 

J . Geophys. Research , 67 , 3049-3062. 

4. For a complete discussion of the cumulant discard approximations 
see Kraichnan, R. H., 1962. The Closure Problem of Turbulence 
Theory, Proceedings of Symposia in Applied Mathematics, Vol. 13, 
Hydrodynamic Instability, American Mathematical Society, 199-225 

5. See Malkus and Veronis, loc . cit., p. 228 for a complete dis- 
cussion of the conservation equations. 

6. See Chandrasekhar, loc. cit., p. 35. 

7. Spiegel, loc. cit., p. 196. 


Markus, loc. cit., p. 200. 

This fact has the consequence tha-- the heat transport, N max 
is asymptotically proportional to R. H. Kraichnan, 

private communication. 

Spiegel, loc. cit. p. 196. 
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Captions 


F igure 1 . Time development of the total heat flux, N ( t) for 
R = 4xl0 3 , IQ 4 , 10 5 and a- 1.5. The system is in the conduction 
state at 7 = 0, with all fluctuating amplitudes except w ]L 
equal to zero. 

Figure 2 . 4.22xl0 _2 w and 4.07© for R = 4xl0 3 and a= 0.8. 

— 3 

F igure 3(a ). Mean Temperature, T(z), for R = 4x10 and a= 0.8. 

3 

Figure 3(b) Mean gradient, 8 (z) , for R = 4x10 and a = 0.8. 

8 (z) is normalized hy the total heat transport, N = 3.92. 

—.9 4 

Figure 4 . 2.05x10 w and 5.16© for R = 10 and a = 1.0. 

Figure 5(a ) Mean Temperature, T(z), for R = 10 4 and 0 . = 1.0. 

Figure 5(b ) Mean gradient, 8 (z) , for R = 10 4 and a= 1.0. 

p (z) is normalized hy the total heat transport, N — 5.82. 

Figure 6 . 4.33x10 w and 9.42 0 for R = 10~* and a- 1.5. 

_ 5 

Figure 7(a ) Mean Temperature, T(z) , for R = 10 and a = 1.5. 

5 

Figure 7 (b ) Mean gradient, 8 (z) , for R = 10 and a= 1.5. 

p ( z ) is normalized by the total heat transport, N — 13.82. 

Figure _ 8 . 8.98x10 4 w and 19.4 0 for R = 10^ and a = 1.5. 

Figure 9 (a ) Mean Temperature, T(z). for R = 10 6 and a= 1.5. 

Figure 9(b ) Mean gradient, 8 (z) , for R = 4xl0 6 and a= 1.5. 

p (z) is normalized by the total heat transport, N = 31.48. 

-8 ' 6 

Figure 10 . 3.22x10 w and 8.57 0 for R = 10 and a = 6.0. 

Figure 11(a ) Mean temperature, T(z), for R = 10 6 and a ^ 6.0. 
Figure 11(b ) Mean gradient, 8 (z) , for R = 10 6 and cc = 6.0. 
p(z) is normalized by the total heat transport, N -= 22.3. 
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Figure 12 . 1.06x10 2 w and 12.93 for R = 10 6 and a = 

Figure 13(a ) Mean temperature, T(z), for R = 10^ and a - 9.0. 

Figure 13(b ) Mean gradient, g(z), for R =_ 10 6 and a - 9.0. 
p(z) is normalized by the total heat transport, N = 5.40. 

Figure ' 14 . Cosine spectrum of the mean temperature gradient for 
R = 10 4 , 10 5 , 10 6 and a = 1.5. 

Figure 15 . The total heat transport N as a function of a f c •; 

R = 4xl0 3 , 10 4 , 10 5 , 5xl0 5 and 10 6 . 

Figure 16 . Maximum total heat transport, N max as a function of 
R 1/3 . 

Figure 17 . Critical Rayleigh number R c for R - 5xl0 3 as a 

function a. The value of a which supports the mean temperature field 

labels the various curves. 

Figure 18 . • Stability diagram for R = 10 4 # a ] is the wave number 
that supports convection, and a a is the wave number at which a 
small perturbation is introduced. The shaded region indicates 
instability . 

Figure 19 . Stability diagram for R = 5xl0 D . is the wave 

number that supports convection, and a 2 is the wave number at 
wl ich a small perturbation is introduced. The shaded regions 
irdicate instability. 
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